{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GPU-accelerated Scalable GP Regression in 1D (w/ KISS-GP)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "For 1D functions, SKI (or KISS-GP) is a great way to scale a GP up to very large datasets (100,000+ data points).\n",
    "Kernel interpolation for scalable structured Gaussian processes (KISS-GP) was introduced in this paper:\n",
    "http://proceedings.mlr.press/v37/wilson15.pdf\n",
    "\n",
    "SKI is asymptotically very fast (nearly linear), very precise (error decays cubically), and easy to use in GPyTorch!\n",
    "As you will see in this tutorial, it's really easy to apply SKI to an existing model. All you have to do is wrap your kernel module with a `GridInterpolationKernel`.\n",
    "\n",
    "This is the same as [the standard KISSGP 1D notebook](./KISSGP_Regression_1D.ipynb), but where everything is super-charged with CUDA. SKI is especially fast on the GPU.\n",
    "\n",
    "The only thing required for GPU acceleration: `train_x.cuda()`, `train_y.cuda()`, and `model.cuda()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Make plots inline\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.linspace(0, 1, 1000).cuda()\n",
    "train_y = torch.sin(train_x * (4 * math.pi) + torch.randn(train_x.size()).cuda() * 0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that. Here we are using the same number of grid points as training points (a ratio of 1). Performance can be sensitive to this parameter (the number of grid points), so you may want to adjust it to your problem on a validation set.\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x,1.0)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.GridInterpolationKernel(\n",
    "            gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel()),\n",
    "            grid_size=grid_size, num_dims=1,\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood).cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/30 - Loss: 1.137\n",
      "Iter 2/30 - Loss: 1.109\n",
      "Iter 3/30 - Loss: 1.080\n",
      "Iter 4/30 - Loss: 1.052\n",
      "Iter 5/30 - Loss: 1.023\n",
      "Iter 6/30 - Loss: 0.994\n",
      "Iter 7/30 - Loss: 0.963\n",
      "Iter 8/30 - Loss: 0.933\n",
      "Iter 9/30 - Loss: 0.900\n",
      "Iter 10/30 - Loss: 0.860\n",
      "Iter 11/30 - Loss: 0.804\n",
      "Iter 12/30 - Loss: 0.721\n",
      "Iter 13/30 - Loss: 0.612\n",
      "Iter 14/30 - Loss: 0.492\n",
      "Iter 15/30 - Loss: 0.382\n",
      "Iter 16/30 - Loss: 0.296\n",
      "Iter 17/30 - Loss: 0.230\n",
      "Iter 18/30 - Loss: 0.172\n",
      "Iter 19/30 - Loss: 0.122\n",
      "Iter 20/30 - Loss: 0.074\n",
      "Iter 21/30 - Loss: 0.027\n",
      "Iter 22/30 - Loss: -0.021\n",
      "Iter 23/30 - Loss: -0.063\n",
      "Iter 24/30 - Loss: -0.108\n",
      "Iter 25/30 - Loss: -0.145\n",
      "Iter 26/30 - Loss: -0.187\n",
      "Iter 27/30 - Loss: -0.224\n",
      "Iter 28/30 - Loss: -0.263\n",
      "Iter 29/30 - Loss: -0.300\n",
      "Iter 30/30 - Loss: -0.331\n",
      "CPU times: user 2.43 s, sys: 100 ms, total: 2.53 s\n",
      "Wall time: 2.52 s\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "def train():\n",
    "    training_iterations = 30\n",
    "    for i in range(training_iterations):\n",
    "        optimizer.zero_grad()\n",
    "        output = model(train_x)\n",
    "        loss = -mll(output, train_y)\n",
    "        loss.backward()\n",
    "        print('Iter %d/%d - Loss: %.3f' % (i + 1, training_iterations, loss.item()))\n",
    "        optimizer.step()\n",
    "        \n",
    "%time train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAADDCAYAAABtec/IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJztnXd4VFX6+D93JpNMIKSHEHpCRwQCI2IHEixYWBEXLKvf3yJs0911FQQVsePq4ir2CJZdRVEEsaICKrIgMpCASAgloQhpECaFlJnMzO+PM5NMuakzSSbM+TxPniS3npl7z3ve8563KHa7HYlEEpxoOroBEomk45ACQCIJYqQAkEiCGCkAJJIgRgoAiSSICfH1AgaDId3x52Sj0Xifr9eTSCTth08agKPz32g0GtcDYwwGwxj/NEsikbQHir/8AAwGwyGj0TjALxeTSCTtgl9sAAaDYR7wB39cSyKRtB/+1AA+BGYbjUaT2v758+dLl0OJpIN46qmnFLXtPhkBnXN+o9G4E8gF5gBPN3T8I4880uQ1i4qK6N69uy/NanMCvY2B3j4I/DYGevug+W1ctGhRg/t8nQKkA7GOv6MRQkAikXQSfBUAGUCKwWCYA2A0Glf53iSJRNJe+DQFcMz3M/zUFkmQUFtbS3l5OeXl5QRqNKrNZqOsrKyjm9Eonm1UFIWwsDB69OhBSEjzurbPjkASSUspKCggKiqKuLg4FEXVNtXhWCwWdDpdRzejUTzbaLfbMZlMFBQU0Lt372ZdQ7oCS9qdmpoaIiMjA7bzd1YURSE6OpqamppmnyMFgKTdsdvtAdH5MzMzyczMbPP7mEwmVq9e3eb3ASEEWjKtkgJAEtDk5+eTnp5OQUFBq6+RmZnJsmXL2LBhA8uWLSM3VyxWRUVFsWpV29uto6OjVe+TmZnJsGHDWL16NatXr2bJkiV1bVOjsX2tRdoAJAHN4sWL2bJlC08++SRLly5t8fkmk4lnnnmGFStW1G27+eabWbFiBbGxsY2c6V9iYmK8tqWmppKcnMy0adPqtk2ZMoUvvvjC69jc3FyWL1/OE0884dd2SQEgCUiio6Oprq6u+z8jI4OMjAz0ej0mk6qzqSqrVq1i0qRJbttiYmLYsGEDY8eOJTMzkw0bNpCVlcWsWbPYsWMHADt27GD69Ols3LiR2NhYkpOTycvLY9WqVSQnJzNkyBDWrVvHihUr+Mtf/sI999wD4HZ8cnIyy5cvZ/To0ezcubPZn9s50m/cuBGASZMmkZWVRV5eHpmZmURFRbFx40asViuTJ08mJSWl2d+HJ3IKIAlIsrOzmTFjBuHh4QCEh4czc+ZM9u3b1+JrlZaWNrgvNTWVtLQ0Ro8ezfLly8nKymLjxo1MnDiRBx98kLFjx9Z1/kmTJhETE8MTTzzBbbfdVneNadOmkZKS4nX8Aw88wPXXX09aWhrJycktanNKSgqxsbHExsayZs0aJk2aRHJyMqmpqV77fEEKAElAkpSURGRkJDU1Nej1+rqVgx49erToOpMmTaob1Z3k5eWRlpbmts05Hbj++uuZNWsWS5YswWw2ExUVRWpqap0WER0d7XbtJUuWMHbs2Lptnse3FJPJREpKCkuWLCEqKorRo0fXbQcxFXDuGzVqlNu+1iCnAJKApaioiNmzZzNr1iyWL1/eKkNgSkoKc+fOZdmyZSQnJ5OVlcWLL75Yt99kMrlNAZwq+8SJE5k8eTLLly+vG32dKrjJZCI6Oprp06fzwAMP1AmFxx9/3O34e+65hzVr1jB69Oi6c1NTU+vunZmZSV5eXt0KQV5eXl3bnPcrLS0lNzeX06dPYzKZyMvLq9tXUlJCbm4ueXl5btdtCX6LBmyK+fPn22UwUPsQ6O07ePAg/fr1C2hHm87oCOTk4MGDDBw4sO7/RYsWNRgNKKcAEkkQIwWARBLESAEgkQQxUgBIJEGMFAASSRAjBYBEEsRIASCRBDFSAEjOajIzM7ngggvcwn5zc3O9tgUr/igNNsfx5wBZGkzSUvT6ML9cp7paPQlGampqnSfgSy+9BIjYAKdffbDjj9Jg641GozM5aHpT50gk7U1UVFSD+3Jzc1m2bBmrV68mMzOz7v833niD3NxcNmzYwJQpU9iwYQMPPPBAO7a6ffB1CpCCSA0OIiV46+MSJUFJdXWNX36aYtq0aSxbtszLH98zgs8z0i4tLY3o6GjS0tJ8CroJVHzNCuyaEXgMsNK35kgkbUNaWho333yzW+Sek6ioKFJSUkhOTmbJkiWMHj2aPn36cPToUUwmk2oyj7MFv0QDOioE7XRUCGqQoqKiJq/VGaRsoLcx0Ntns9mwWq3tcq+srCxef/11+vTpw+jRo+nduzfbt28nMzOT7du38/DDD5ORkcGkSZPo378/ffv25eDBg5w8eZKDBw/y6aefkpuby/79+8nNzWX79u11IbodTUPfoc1ma1ZfAz9FAxoMhnlGo7HBkmAgowHbk0Bvn4wG9A8BEQ1oMBjmODu/NAJKJJ0Lf6wC/NNgMBwyGAyn/dQmiUTSTvhqBFwPnL0WEonkLEd6AkokQYwUABJJECMFgETSCtqz3FdbIrMCSzqcF7495NP5d00c0Oj+zMxMduzYUZebf+PGjc2usON0DMrKyqor/gH15b5cq/p0RqQAkJzVqJUGa+7IbTKZKCkpIS0tTbWM2NngISgFgOSsZtWqVV7uv/fccw+5ubls3LiR5ORkSktLiYqKYsmSJdxzzz1s3LiRhx9+mB07dpCXl8eGDRt48MEH2bp1KyaTyavcl/NazpJgJSUlbtd64okn6mr7OdsyevRot3M6KjJR2gAkQYezjNcdd9xBWloaq1atUg36cQYJpaWlMWbMGADVcl+eAUVq11qyZAmzZs1i2rRpTJo0yeucjkJqAJKzmunTp/PnP//ZbduGDRsA6ir8mEwmn4N+XAOKQH164JxGOCsJeZ7TEUgBIDmriY6OdisNVlpayujRo3n88cdZtWoVsbGxTJs2jby8PPLy8upKbWVlZVFWVlZXCmznzp1kZmaqlvvyLAnmeS3nec8880zdqO95jmvNwfZElgZrBYHexkBvnwwG8g8BEQzUmai12igsq6aiupb2EnwSSSBz1k8BaixW8k5VknuykmMllZitNgBCNBoiw0Po3i0MQ79oYrqEdnBLgw+b3U6NxUZ1rRWbDRQFNIqCRoEwnYawEG1HN/Gs56wWADmFFXybU4zF0eldqbXZKDljpuSMmf2FFYzoFcm4fjGEh8qXrq2x2aGsuhYrtgY1sSqLFZ3WSkRYCKEhQaWo+oTdbkdRVLV9Vc5KAVBrtfHDwVPsOVHWrONtdju7fy0lp6CC8ckxjOzdcBJJiW/sKyjnZKUVbWkpXbtFNfqyWqw2TleaCQ3REKnXodU0/8UORux2OyaTibCw5mdaPusEgKnSwrpfCimuaDpRpCc1tVa+P3ASU5WFSwbGtUiSShrHarOz6cBJ9pwoQ6dosdSWEGUy0dxvWKMohIdqaS8ZYLPZ0GgCW/PwbKOiKISFhdGjR49mX+OsEgCnK82s2nmCaotv+eZ2/VpKlcVK+tDuctTxAxU1tXzxcyGF5dUAWOwK2/Jr6dq1a4uuE67Tct2oJLp3808tgcYI9JUU8E8bA1vEtYAqs5VPdxf43Pmd7C+s4NOfCzDXetsPJM2n2mJlbVZ+Xef3hSqLlY+z8jluqvJDyyRwlgiAWquNz34uoLTK4tfrHiupZO2ufGpVjIiSpnE+l5JKs9+uWVNr5ZNdBeSX+i5QJH4SAI604B2C3W7n6+xiCsra5oUoKKtm/b5i6TfQQsRzKWqTjlprs/HlnkLO1NT6/drBhj9qA6YDrwGNB2W3Ef87VMKh4oq6/8tKinln8b1M/eMC1r66mPSb/shbj92F3SY6sN1ux2qpQRMSiq3WjFYXhqIoju1mEnonEx4Rwe0LlxIZmwDAgaIKorvoGJ/sHRIqUWfTgVMcKj7jtb2spJi3Hv0rdpuN6X97hPeeuY/CIwfdjtHqwtBoNMQl9UYbEoI2RMftC5cC8M7ie7l1wb8gNoF1e4u4flQSGmmnaTU+awCOxKC5fmhLizl8qpLMYybKSop5ee7tlJUU88Ubz5K3x8hzd04nb4+RZQvnYKmuotZcTa25GqtFrA7YaoVaarXUuGy3U/xrLkf37Wbtq4vd7rX98Gn2FZS390fslOz6tZTdx0spKylm6d9v4oW7b6KspJiykmKe/fM0juXs5tcDe3jh7pu8Oj+IZ2KpqaLg8AGOH8zm6L7dPH5bGmtfXSye7V03UlZSzAlTFZsPneqAT3j20GlXAarMVjbsKwZg/YpXyNtj5LFbJrgc4Rjxbc75exQwBRgLnAAOOn5yAW81dfcPXzH3qnMAhYXvfktkbAIbc4qJCtchXYUapqi8hs0HRaf84o1nOZazG8Dj2QisFqdtYBgQCZwETgGlOJ+fE7vVyu4fvgKgvKSYx26ZgKLREqoP592PPuOaiePb4uOc9bSrAPBnabD1B0z87UoDtZbGDEzdgN8BvwEmAGrBHRbgDWARUKiy387KZx/k5gXPAvDBj4eY2CewRUBHlQazWG2s/aWE+64f38RzAUgCbgZuA0Z67KsFfkQ8k40NXsFus1JTWcGf5/w/Bq/fQE3ZKe68805eeukln5fHAr28Gvinje0qAJr7UJo6bm9+GafMWha89TVrXnqMPVs2qBw1BviQ+oLFVuBbx088MAgY6Nj/B+AW4J/As0Cl25X279jMw9PHEaILZfEnmfxcUsvI4YG9RtwRa9jr9xVh1Ta1Rt8TeBm4Bup0qZNAHhDn+IkCLgY2OH4eALY1eMWiY3mMHCKes0ajISMjg6VLl7b+gzgIdD8A8L2NnW4ZsLTKwqYDQsVc/H+XN9D5/wBsQXTuTMQo0x2YBDwG/A0xHRgMnAOsBSIc+/YDaar3rrWYOZG7j2OmGnb9WurHT9X52V9YwbZfcln695tI7DeAmMReKkdNAHYCUxECebXj757AOIQdORqhuS0ATiOexY/Ae0B4k+2w2WxkZGSg1+s7LMa+M+GP2oDTxS/DdD+0p1Hsdjsf/LCH5//xO+ZfO1pFzewKvAu8CoQhRpoLgP8CJQ1cNYf6KYIR6AV8AcxUPfqFu29m2YLf8+X2fZxshbvx2UhZlYXv9p9k/YpXOJazm+MHszldeNzjqHnAeiDR8bsvcAPwCWIa5koF8BSQDDzu+H8msA5hK2iamTNnsm/fvlZ+ouDB5ymA0WhcBazyQ1ua5OfjZax47Tny9hjRhHjO58OAL4FLEC/MbOD9Flz9e8Qo9DRwL2LESQSedzuq1lzDrwf28PBNl9Fn0wF+O7YXOm2nU6T8ht1up1diPBZzQ8JQjxDKzvTZTwAPAd7OVbowPZYaV4NsKbDQcf43wKUIm8CViGmDN4pGA3Y7kZGRLfKJD1Y6zZsbFRXFuIGJbP18JQC2WtdRQwHeRnT+Y4CB5nT+UH048T37umyxA3MRAgDgOWCx13lOZl06iNggVzN/yS9n/ptfEauq8muo7/wm4DrgQdQ6v6IoPPnxDha++x1xSX099u5D2AQOIFZxfgB6q7Zn6HmXEJfQnSNHjrTyEwUXnUYAvLDmByHdVVkMzECMGFMQan3DDEq9gMR+AwnrEsHA0WL5yP3aSxCrBxZgPsIwqM7IS6/i2Ong9E2vqK5l/NA+PHbLBEq8VH6ApYjOfxrRgT9VvU58z748+M63AETGJmCzqcVzHEEI+F3AUIQQiPc6Knvb95wsKiSpl7qAkLgT8AIgPz+fceMv4vF/3MG5F6arHPEH4D5EZ70B2KN6ndjEXtz90kdccPVMwsK7cu+ra3no3e+oMJVwwdUz+fsLH3oIgXcQI1YNcDfwR9Xr7tiwlkFJMUQFoSbw3YGTLHjzqwb2LgD+gvCxmAr84nWEM9zaZrXWeV0C9BowjAuunsmAUecTHuE65y8ELgN+AvojtDz1Jdm33lguDYHNIOAdgR57/Al2Z+0A4Ci7PfZeCbzk+HsOYsmoHkWjZfRlV3LNHXPrXrBpdy50O+b2hfVz/Af/u5Glf5tB6UmnP8A6hC3hP8ALwCHEXLQeXZieERemcd9Dj7XuA3ZS9hdWMPPCwQ2s998OPIlQ9W9BjNbeaLRaxl0xnfLT7vN512fy9mN/o1tMPFs/d07pShFG252IFYInEFqaOyGhYUR268a6L79o4ScLLgJWAERHR1Nd3VggSRJilNYCjwJveexXsNus6LtEuI0ujREZm8CwcZfx4xcfuGz9LzAEsRb9IWJVIbtur6Wmmsxvv+DXec+QX1pNUpS+WffqzOQd/ZVpN8xsIEDqYmCZ4++/Ipb6vAnRhbLgra+bfDZOYVBSeJwco1OQ5AO/RRgE7wO2Ax+5nVdrrqHkVA2vZbzOSy++0JyPFZQE7BQgOzu7iYwsbyCcRtYhPMbcGTBqHBdcPZPy0y3zFXdOCbRuqwwLEZ0/CvgM77mnnV8PZXPF5ZM5fvxEi+7XGfnHAw9z8Gcj1lrP5bsYYAViXPkXTu1Mpw8nVN8FECOzoiiMnnhtswUzwB2PvUpUvKtV/wfqjbVvIuwC3ixf9rqcCjRCwGoAMfHdOS/tWrZ9s1Zl758R6v8p4Peq5x/atY0jezNZ/Elmi+7rHHHSb/6ji/+6HaHW9gfOQywRXo6rv/pzd94AwN33P8wHb2e06J6dhaa1steBPsBWhA1AYKmuAkXhgqtnMn7Kjfz4xYecLi5o8f37DB5B1ZkyzFVOT83ngfOBm4A1iBWCStVzN23a1OL7BQMBKwC+23+SM2e8w0mFOv4vx99zEOqgO855+TV3zAVECvC+seGkxHelZ7SeMzVWSqssmKosFFfUcOSU90sTGZtAVHwPSk86X9QqhDErC0gH7nFpRz2frPwP+pX/Qa/Xdwp/8pawZ88vDBo8yCXAypU5CCNsKcLH3z1Wf+yka+vsL9PuXFj3bCP1Oob0iKBLqJZqi40qi5Vqi5WjJVVe2Z32bd+kYnO4AxFLcA7Caegfqm1ftmyZX9yDzzYCUgAcKKrgQFEF+7Z7Su0QxLw/HLHuXz+/1IXqsZirURQNteYa9F0i6NEjiUsHxdE/roubs05UuI6e0fVz9fzSajYfPOWVVKTslGfwUj7wfwhPwScRcQU7vNqvKBqys7O9tnd2jlvCOffCdHZv/tpjz3CEzwSI1ZLDXufu2PAJOzZ8QogulH99voshCV04f0gvt+fgirnWRuYxE1nHSutqOSx462tW/HMeh3b/5HJkJcLVexvCxftDhAbiTkZGRp2L8NkmmH0h4GwAlWYra7fuZd6UESrSfiHCyecwwsBUj8VcLZbzXvyQ8VNmYKkwMcPQi0HdI5r01EuK0nPj2F5cdU4iUeH1c/8H31GLRPsSoXrqEFMB78SWYyZdQ7G1S1MftVNRcsZM2qhklc6vRyzHhSPsMp4OWGKpTxemJ3Xi1Tz9wXfMMPTiouTIBjs/QGiIhvOTY7ltfF9S+0SjKAqRsQlUmNRsOjsRHpwaRxv0dfd0ogvT89sZM6R7sAcBJwC2HCnn8/+8hN1u9wgoORcxr7QhJH59zn+NRssQwyVMu3MhvVKG8cTTS9j05cdEhLVMwRnYPYKZht70jglnwXWpqjHsgvsQDimDEM4u7tRUVbI1t4SK6no1OD8/n/T0dAoKWj73DQT6JSVQq+ru+zDi2eTgKZQ1ISGAnZDQMGrNNfSIj2H25aNbVIUpPFTLxQPjeGBqKnOvOofCo4caOPJRYC/CGCiMwq5uxZaaaj5YubLZ9w0WAkoAREdH8/u0kXXuvvUBJRrE0pIOEeDjvq5ss1nJMf7AgqljuGpEIhcOiGt1mqjQEA3XntuD1z7dTOqEqxvwPqxBGJ6qEEbI34rNigZFo+WCa2byr7tmcOHFl9R1+MWLF7NlyxaefPLJVrWrI9lXUM5fnl1B16hYj+9jLMISb0UIZXebja1WCEC7zca0m28n1Fze6riJffv2ce310+tGdUWjQacPJ7p7kuOIGsSzsCLcuQ1u5zvbPe788Z1SCLfVABJQAiA7O5uwLmq54u9EBOr8CtzvtVdRNKROvJrPfzAyIKFluebVCNFquPmyc+mZEAN2O4oiviaN1lWjyEZ4CIIQSt3BbsNus/L6/XdwLGc3+/dkkjJgAHq9noyMjE4ZqlpltvLDwVNs+/IDzpSWuBgAdcByhB/GcwjvPG/GTLyGDdt28e7yV1npwwiclJREj/gYai1mQkLDsNtsWKqrMBW5GoG3Af92tOlNXBPAONtdVFhA//79O83378RzAPFXJGrACAC9Xk9ycjI1lZ6W/74Iby8Qy39qefns9EmM47JRg/zWHkVRCDGXc93M29BohbupzeqZhfY14CuEP8KLqtexWd0t2eHh4Z0mVDU/P5/YqAjuSh9ap5XVMw8YhfCOfEj1fEVRSE6KZ/w5Kar7W0pRURFzZs/my282cv7lvyEqPlHlqIcQQUMjaGhFAKC6urpTCIHo6GjVAaR/UvN9KBojYATAtm3biIiIUNnzCiJZx4d4BpNoQ3R0iYzhyum3olT5P0HHypUrWfnma7z55RZSJ1zdwFFzEELpRupDXt1J6NETRVHQ6/XU1NR0mlDV+xc9puLtpyBy+DldqmfT0Nq73W5nzfv/8Vt7Vq5cyfPPP88l54/lnbffZPj5E0SLFNfXuAoxUOBoo2dkoUCj0XQKIZydnc2MGTMIDxfJUJzOccPHT/TL9QNCAERHR3P++edTUVHhsWcmIrrvNHCX13nWWguW6ko+fOs1n9TLprjhohHExoiCoe4vG8BRhFEQhOebd+rw4oITKIrCpk2bmD17NoWFarkHAwfnqPPef95Q2asg7DFhCMefb1WvodFomDlzJrm5bZMwuk9MODpzGRdcPZMRF3kGia0HViJWaJ7zPhmROWjoUHXvwUAiKSmJyMhIqqpExKnNMZXZ/cNX9O3b12ctJiAEgFPKuRNHvYV9Lp4JO3Vhes5Lu5aff9lLSBsn5NBqFMJrK5gw9WaVlw1EBqLvgR6IOag3NpuNSy+9lOeff75NhZU/yM7O5uK0K0G1dOefgAsRmZXnqp6vCxV5Adta0/n60zUYv1nDz15LkyDU/3LgeuAqr71R8T3YbPQMLgtMioqK0GrVox59LVgTEAIgKSmJVas8kwr9G0gAvkOs7bpjqalmx7ef07+PWiIK/7Pqww94781XUew2Rk6YwhDDJS577QiPtCqENdz7hQMYMOScTmGB1nWLxaqPwjM1t8jd51zFuBPh9aeC3dZums4ve/eSeslkldWaE4glShCRnO4+B8PPv4xfTBpstsCv+LRy5UpycvbTvVc/t+39+/cnJ6fx3BdNERACACA9Pd2lHPdViIQclYiOpf6QJk+e3D6NcxDTJZRXXlxKaVE+EVExKIqCVudc0z6IyHYDQiNwX43Q6kLZu3snix4N7LBhm81Ocs/uKkY/EOp0JCKP3xrV86+9/kYOHDjQbppO3969GDGgTwPuyUsR+SEGUD9NE2z9fCUHj/zKhZdN7BRC+cCZUCyOZVUc/SQ/39sNvqX4JSmowWBINxgM83y5ztq1ax0CoBvCug7CiKPu+KHRaFi7Vi1QqG15+5XnOJqdRe6eHaAoLsUtQHQQI8Lw9KjbeVaLGbvdztsBnqhi+5HTzFdN8nEVwtB5BjV7jJOkhJh2N3CWnT5Fz959GXnJFYy85ApCQp2pyWsRUxYQOQPqq9eNvPhynrvrRnYZtzFv/oKAdtI6WFTBnhNl9BowjMR+A8BuJzwikpqaGp/9SnyKBXAWBTUajesNBkOKwWAYYzQad7bmWtHR0Q4Dx1OIiLKf8EzI6SSiWzcuuvDCVra6dXhGwjmdlBRFAUVxjEA2xKrAdoRf+rsIN1V3lAC1QBeW1ZA2KlnF4y+c+sQrixCGT3ceffNTjm39tEMMnCtXrhQrDi6lw997Zj47N34KbEbEjdyOeJ+uAXBzaf7g/fcAGDBgQAMBaB1HaZWFjTknWXBdqptrfFWF8IT1NcbBVw1gBiLbI4gaW2oWsmaRnZ1Nz543IZZwzMAshFeXN/Fxce0++nsux4SHhzN20jWMmXSth/qZidAEtEAGaimrxk66lrDIwCo0Wmu1sT67iAVvfuWRKBXE2noywv25Xig7HbAefX8Ts6dexgtLl3aYgVNRFCYNSSDEYQswV1cRGu6Mx5iHsFdcDVzb4DWsVmtAaWdWm52vfimiuKiApJShjLgwzUW7EWi1Wp/8SnyNBozGPeF+XGMHN1QabNCgQdTUKFCX8utJ1HL7dY2K4frrruFkcXGzyoz5E61Wi06no7q6mtDQUKqrq+kZ240DR/MZOWEKZScLObzHGRm4CJiOcJW9E9dOExrelYqyUlb/dIirh8WiUfxf2balI0FhYSG33/Encn7ZpZLk4xxE6LMNkX+x3hnKbreR9d2XPPX0EiwVJoo8V3H92MbmMjQGth09w/R/PMkjNzrrBRYhppNLEc/iG9TqQYKYWm7evDkgIgZ/OlpObsEZvnz7BX7d/zNVvfq7aGcKYMdqtfL++++zZs0aDhw40OJ7BERpsH379nH33Y/x8cd7ED7d3qm4NVotF190Ectef71tG9kI5eXlzJkzh+uvv541a9ZQUFDAU6+8zYGiCuZdfa7LkWcQCTE/Q8Sor0akKwdz1Rn2/fQdf71mHBmpY1n74fttMmduScmoBQsfZu8uo9foIl6yDIRL7Suol+eyc+mIfi4G3LZpY3NJSLBz0nKCgrJqHnxnI5+9/gx7tm7AUvMyQqschbAHPKx6fnx8PPHx8Wg0mg4tDZZ78gwzJw50U/tP/prncsQ7iPD0+UAtOTk5rWqvr1MAE/WeL9GIFD0tJikpiYQELYpyA6L4g3eiyejYONaubpf6Iw3i9EQbPnx4nZV74uB4uoaGMHjMRR5Hf47wXoxAzU3YWmth1/Yfufe+BV772gunw8+7b4tlVu+5/x+oX/P3bmd8z35k7d3fqs7fViiKQtrQBLQaET6ctelLR1SgFaGNgVgR8HZPDusaQVFRUYcHbJ2uNPPR5l8I66rmGQvCCHszwt7Ul1vUgKMIAAAX1klEQVRvvbXVg4ivAmAl9d9kCsIFq1UUFRUxbNgw4DRh4V2Iju9Br4HDGZs2ldEXT+Zw3mEfm9o2hOm0/H3KSJeEla78DTH3vA6RLcebVSvf67B55969exuptdATYZAF8cJ5r/mfPHGEIcl92qh1rSe2ayiGfjEAHoJ5MyLDsx41D8GaM2IOk5GR4Rcvu9ZgrrXRt0cCC2deyhmTWjm7CxF1KwD+H5DLO++80+q2+iQAnBZ/g8GQDph8WQFYu3Yte/fuBUQ8velkAccP7iVUr+e9998nNCRgXBa8yNm3j0uunKrSmfKpX39+AaEkqeOrR1dzcQ0rPVKtb6DWAoj5chSicKp3Zt9QfTgTJkxou4b6yJg+Udw/dYyKYHYaBK9FpHjzJjRMz3VTp3bISk1cbEwjZda6Ax8gpmRLcGZCnupDW33uVUajMcNoNK43Go2tzoTpaWF3ZevnKxnWOy5gLLNqJCUlseWbzxpwRslAjDxJiPLj3sT17Mvqb41t2MJ6nGGl9y18hB1HTCoZfqBeYymnXm12Rx8Wxrp169qwpb4RotWw8cdMlZwOhdQ7bL2I8DtxOU8XisVcQ7ktjMREtWjDtmPnURPz3/zK0WbP1SMtIttSL2ATrrUQPvvssw6bAvgFpyuwM+DBk84QPpuenk7f/imE6Dyz3dgRczWz4/clXufarFayy7RtWmLMM6x05X/f5N6rhqsc2Y36Nf8HEDkYvCkrNQXUkpka5w1PoXtcNHhpVy8jDJq9qQ81F9Q6HLbyC4vYdLBVJq1WsX1vLrdOEz4KYV26Yvcqj/YkMBGhVc7AdTUmPb3Vq++BIQBAfIiBAwd6dCBF5ILrBOGza9eu5YrJaQ1Uysmm3oc+AxFJV8/pwuOUnipi3S+FlFV5LsP5B08tKyQ0jJjuPVUs/08hOsZP1AsCb7RarU+qZ3vxzZr3VKZXNkQYswWxWjO+bk9UXCJ3v/QRV8+exx9v+g3f7Njf5m08bqpi/kOPkvfLDtaveIVtX3oau69BTF1qEdmn3D0Wv//++1bfO2AEwNq1a8nLy/PoQHbsdjvLli1r8LxAoqioiJk33UJ0QpLK3sUIQTAUtaxG61e8QrXFyud7CjHXqk0lfCMpKYkuXSOorq4mJDQMq8VMqD7cw5V5MvWOWHfgWsU3PMJdVbZarSQmJga8YD506FADBWZ+RqR111C/1AmlpwqJiI7j+w+Xc/iXHTz55JMcaImDQwuJiopiQI8Y/vfZ+9jtdrZ+vhK73fX596a+6tUCxHTSHV/sRwEjAEBoAXFJfeq0AKeXU1vFlPubr776ivffexdTsVqQhnMKAOJBjnHbu/Xzlcy96hxmTxrOx1n5VJnVvSBbS6XZyq4DRxk/ZQZ3/XsF46fMoPDoIZeXJ4r6qMtFiA4iuODqmYSE6unesw833HAD06dPp3///gGf1wCE4Jsxc2YDex9FBHGdS32VIXjslgkYv/7I0SHf59y+8URFRfm9bacqzCz8zzekTrjaLYNxPSGIeX8cYll5idcRAwcO9CkiMKDqAqxdu5bLpt7CT1+vRq/XYzabO4X67yQ7O5v58+fzwQcfNCCVNyOs639FLEeNRTg+ifiAwWMuoqqilANHjrHaamPq6KQWZzZWo8ps5eOsfG66v37pa/vXnpb95xGjzVbgmbqtEdFxTLtzITP/tohbzu/jl/a0N5VnzjBo8FAO7PecrlQjfB02INydPwLcVX5nkZmpc+axr6CcoT264Q9OV5pZuyufsMg4wrp0dctgXM8jwEXAcUQsg/c7VVtb61P/CCgNAKCq/DSz7pjdabLnuOLM3qIoChovK66T+YgU2ucA9aHBdpuNHOMPHM35mfUrXqGk0syqnScwVfpmEzhVYWZ11glOnalfWiorKUYX7lq3YCriBat0/Bbah6JouPulj3h57u0M6lbbKTs/CAeu4cOGcOWNv1PZuxGRQFSPCN7Sue11FpnpGh3PN9lF/O/QKZ+XbHMKK/jAeJwzZmHIqzCVEN+rv8dRlyOmilZEBmp1g+TRo96BWS0h4J7osowMRg7sDcDzz6tHAwYyRUVFzJ49m5ycHLZtN6LTd6HsVDH10tuZNGQLwsf+E9zmdY554NbPVxKiC8W67mcmD0ugb2zLCo3YbHZ2HjOxLe80NpcX1jOqTBQ6dYZfz0ck1BQMHnsR61e8wuFfdrDmjaVMeqHzVtn96quvGqlr+HdgAiKV+KO4ez0qbgVmdx41cejIcd5+4h+8++47LRp9LVYbmw6cYm9+fU0L7+cBYqnvv46/F9FgeXWNxufpccBpAD0im180IhBxuguvW7eOoqIinlm1mYXvfktckmuE3U8Io6AGYeBRT2W+4K2vqTTXsnZXPh/vyqe4vHmpoEura/ko8wRbc0ua6PyK4/6JiNx+7i7LOcYfHEYpO6+/3rmr7DpXQcL03r4mosjMrYjRdh5CGAjsdht7trg7uL750hK2bPkf/7j/4WYbbAtKq/lgx3G3zg/iGbsnnA0FViGcfr5BLS7GybRp03yeHgecBnA2odNqmDA4nnMuVytz9igiPDUVeBYxF61n2PkTeGfxvdy64F9ExiZwrKSSlSVVDO0RwfCkSCLDQ+gaqq3zwzdVWjhSUsnhU5Xk/HqK8C7eGsOoS69kx4ZPXLYsdLThFKLmYcOqbXh4OFOnTuWpp55q8JhAxjk9s5hrUBQNdruNUH045mqn78UWxJTsYcToOwoooVtsAjPvXczLc2/nWM5ut+e4esVbrF7xFqGhYRwvOkVEmNYtLsJUaWF/UQU5BeWYVJZ3y0qKeWfxveTtcXUCex6xLHkEofqrC5ie/QepJNFtOVIAtDH94rrw+meb+etNUygvOemyx4JIe7YDsTqwBZG4QpC97TtALA9Ou1Pk3bdjJ7ugnOwCURtBq1GI1Ouw2+1uL1hpSTFvLnqoTnioq5lXIdRLGyKwxHsuqe8SQU3VGcLCwjpVOvOGOHr0KPHx8cTExhPRexC5PxsxVx93OeJxxFLoRYilwemUlxTz+v13ABAeGU2ovgvmmipqzTVuVajf2noEBQWdViE0RINGUSirVrffODt+bGJPj87/e0Rx1WpEivmGHZEKjh7ix00bfPg2BAE3BTgbmXrBOYxU9bn/BeGIAiKPoMHrCOfy4ILrUr32WW12TleavUYX5xr2F288y8tzb6/Ln1+f5TcZYfDSIKzfau7AUF1ZwZw5czqlQVaNfv36UVRURM6+veRmbaV7nxQPN2ErYipQinCF/rvb+VVlJirLTdSaa+rqHeq7RBAZK4p02LFjttqoqKltsPMDPP67SeTtMXpoYwbqHa/+iFomKRAOXFNvuNFvS+NSA2gHuoRqCbWUEZPYi5rqSipLT7vsXY7wCfgzItHmWEQCi3oURWHBW+qd1BXPkd79BQOh4usRy10xCANkw6GvV1xxRZ0htjMaZJ14pnMDOH2yiNMni9CG6FBCdC6h0IcRXoIfINbdj6IWDOU8ftcPX5F+85/qhEBjqGtiICIvP0I8m5dw1QTV7psYJ/Iu+iMpjtQA2okv167m5bU/8Mj7m4nr5Zly6+8IS29vhAHIfSlqzKRrm/eCeRmUPAlBzG9TEQ4wt9HQvH/gwEEdknS1LcjOzua6665T3WettWC1eI7WHyJWAjSIxBsXqJ4b36s/VeWlrF/xSrPaob58GI+Iou8L/I/6epPq3HrrrX7VxKQAaEcmDhaJKsxVVYRHRLrssSAy7v6KCBZyX24TyS2bJjI2oYHiqiCiyf6LSFVWiphjNlxO7eDBA53a6u9KUlKSamSfVqvlwiumMnjsRYxNm+ooZ+7kKcS0LByhKXnXnTx5/HCd+25D0zRX7nrufRS3ArNRiNqSwxCel9ch3gV1brnlVpYtW+bXvItSALQj0V10jOsfw0PvfseAkeMIdVuSKkRUsXF6p71K3eNRNJSVFLP07zfxwt03cfzQPl6eeztlJcWUlRTX/Q2oBJLguM7biFJrZQgnk59VjqunMxUxbQ5FRUX079+f5OTkum1Wq5WBPeOZ8/hrzLz3ScZd7lnb8U6EC2488CWi8pM6cT37NjpNKysp5rk7b8BeV2C2i+PaYxDeh5NxT6/pzpChw6ioUCuM6xvSBtDOjOkTzcRz+zeQ9MGIGJk/QgiBKOA27DYLX7zxLMdyRNLUl+65hVpzTZ3qedgRRTbtzod48J2NfPD8w+T89J3jmhqEp9stiPj+K2molLeTzlbEtDk4R80ZM2Zw+eWXM2vWLJYvX05BQQGGftH8dPg0FV4ZeKyI0NvvEbaZnxBek5le1z914iiP3TKBEF0oiz9x3+89949A2HsuQiz3peNZ+i4svAtDDJegAKeO7GPokMFtknFZCoB2RqNR+CnzZ26b/Sf2/PitShKRL4ErEAlFZyIq8Ux3M+g5/cZdq/c4vQfdSUAU8JwKVCA6/1aVNmnpFpvAuSNGcM6QAW6d42zDtRM5DZtWm52DRWe4feHzKp31DGLJ1NlhNyN8Jj70unZ4RCR//Oebbtu8rzcGEeAzCBHWm44zYWwdisKg1Av53f3Pcl7/GMYnt10KeSkAOoBhA/oysF8vft7SkBfZD4jkD+sQ1ZG/RVS48R55XOkSGU1lmTOd9TTENCIBMde/BuFr4I3NZuWCCel8/E592HVntvq3FK1GIX1Ydwb07t6Alb4YmITIjPx7xArBowinoXrDXlVFGf/+yw0kjzCQftMfef2BOzyucxciBDkUyELYfQ66HREa3pX7ln1OZGwCvaLDGefIbdhW+MUG4KwQJGk+tspSYhN7Ep3Qo4HEnDsRBsFjwPmO//9LQ/XuAUfn74GwXH+E6PzfIMJdvePIXfl69YqWf4iziMTIMFZt+KmRVRRnsZq/I6YGDyHyO8zGs/Bo3h6jR+e/EGFIXIro/C8ivP3cOz+ItPGL/+9ywnVaLh/eHY2mbTMu+6M2YDpq+pCkUVauXEn2vhyGj7tMJWWVkxxgJGLUqEE4qeQgvNT+ikhsOQIYhxiNfkKkjLoFobr+GTGd8FAxXVA0Wq65fnqnybnQllxhGExcTFNx/88jpgSHgSGIZ3EE4VX5G4SmcJ7j53FEwaz/IZ6VCaGZ3YUzDNyTqPhE7n/7G648J7Fdoi99voOjLqB8e1pBVLiOsNpyxk+ZwYiL0lnz4qOcPOHpkmsC5iIcRJ5AuO3ObuSqVYipw1waKqxaj4LdZqVnQuxZY+zzBUVRCLNUcPG1N7Hjuy+pqaxwFOGspNbNkegbYCBiSXUuwkD4cCNXPga8hxj5GxHGisLw8ydw1XlD6B2jFrTkf9rVBtAcz6VAKMnUFP5s49uvv8p3h0rJPVVFvxEGTuUfQxuiU5mLHkaM7M8gDEfJLj8K4qX8AmEvUA97DQ3vQu/B53Js3y5qLWYGjxzLuHOHcvTo0XYvtRaoz/mVl1/kcEk1X+z5E3q9UO3ff3oeB3f9SK1b0g4roizGSuAyhHYWjzDadkMs821BuFxvQs3hStFo6RoZjbm6ip4DhxPfqx/mspP0CjO3W18JiNJgrT2uI/FnG6+Pi+ejnSeoqShj/JQZjJ9yIz9+8SFbP39f5egsx0/ziU3s5bVG3U2vY8bYXoSHNpS4pO0J1OfcvTucKDNz7IyYf8965CUevWUC5eaaBqZr3zt+WkaITsei9zbV/R/XNYwbx/ZEp23+zNzX77BJAWAwGOaobM41Go2trgIkcUen1TBlRCLlj75ItUVk45l250L2bN3AGVMJmpAQlbJdjaMN0Yly63Y7PQcMc9+nUbjqnMQO7fyBzvh+3VAKbBwtqQTgoXe/Y8F1qejC9NisVmqqKxux3TROiC6UyNgEt+cSFqJlyojEFnV+f9CkAPCl4Iek+USG67jynETW7sqv8xl/6N3vWnWtM2fO0LWrukuwgsLEwQkkRnqmA5e4olGEkFydeYLiCiF8PR18/IVOq+HakT2I7qJr+mA/449VgOnil2G6H9oT1PSJCefiAY1WWPeZSwfHMSzJP4ktz3ZCQzRcM7IH3drQGq/VKFx9bg+SotSyArc9/lgFWIUIYZP4gdF9oqiutbL98OmmD24hlw6KZ2Qv/6e3PpuJCAvhulFJrNp5gppa/6Zqd2oZfdrJ4q/ahg67s6RBxifHcvFA/2oClwyMZ1Rv2flbQ2zXUKaOSqJrqP80AUVRuHx4d5LjG4rebB+kAAhQUvtEkzY0wS3HXGtQFIVLBsYzuo/s/L6QGBnGbw29SOzmu6qu12m55tweDOoe4YeW+YaMBQhghidFEqrV8HV2EVZbyy3OMV1CmTwsgcTIjplfnm1EhIUwLTWJjTknySlsXWhur+hwLh/ePWBqLARGKyQNMrB7BImRen7MKyGnoAJ7I5l7nSgojOodxQUpse2+rHS2E6LVcPnw7vSIDMN4xFRX3KMpFEXhvH7RjOsf47NW50+kAOgEdNOHMHlYd0b1jmLLoVMNlhEP12npF9eFRF0YIwfGt3Mrg4uRvaM4p2ck+wsr2HnMRMkZtShCiNTrGJwYwZDECGK7Bl7NCykAOhHdu4Xxm9E9KauyUFpdS3m1hdKqWjQK9IvtQmJkGIqitLtbb7Ci1SgMS+rGsKRu/Hq6ioqaWsy1NsxWG1abnb6xXTpsea+5SAHQCYkM1xEZrkPkq5MEAu0VvONv5ARRIglipACQSIIYKQAkkiBGCgCJJIiRAkAiCWKkAJBIghgpACSSIEYKAIkkiJECQCIJYqQAkEiCGCkAJJIgRgoAiSSI8TkYyCVt+ACj0Xifr9eTSCTth08agKMu4HpH6vAUx/8SiaST4OsUIAVRpwpEFcQUH68nkUjaEZ+mAB5FQ8YgCqU1iKwN2D4Eevsg8NsY6O2DAKoNaDAYxgA7jUbjzsaOk7UB249Abx8EfhsDvX0QOLUB06UBUCLpfPhcG9BgMMwxGo1PO/5Ol0VDJZLOgz9WAf5pMBgOGQwG/9eykkgkbYqvRsD1QIyf2iKRSNoZ6QkokQQxUgBIJEGMFAASSRAjBYBEEsRIASCRBDFSAEgkQYwUABJJECMFgEQSxEgBIJEEMVIASCRBjBQAEkkQIwWARBLESAEgkQQxUgBIJEGMFAASSRAjBYBEEsRIASCRBDFSAEgkQYw/SoM5C4NMlpmBJZLOhT+Sgt7oyA04xlEfQCKRdBL8kRTUmQY8panCIBKJJLDwV2WgecAfmjpu0aJF/ridRCLxE4rdbvfLhQwGw4fAbKPRGPhF1SQSCeBjaTDnnN+h+ucCc4Cn/dtEiUTSVvhaGiwdcM77o4Ht/miURCJpH3yaAhgMhmjgt45/xxqNxibtABKJJHDwmw1A0jEYDIbpgAkY4yzS2sBx8xrbLwl8DAbDmIZW2pr7Hnjil1WA1tJUo1v7odqxfU77yICOcIJyscGsNxgMKQ29IA5/jcl0gH2mGd/hGCAFwGg0rmrn5jnb0Nz3MKWpatltheMZvgYMUNnXrPdAjQ5zBXZtNGDydCJqan8AtC8dWO94IVJcPCLbkxmIFxOEEbYj2tAgzXyGCxwdP6UjHMma+R7mOvbndpSzm/P+Dexu9XvQkbEATTW6o1/upu6f4rIt1/F/exMNlLj8H+d5gGM0WO+5vZ1o9Dt0jKzbAYxG49Md5EjWnPfsn47fgers1uR70BAdKQCaanSrP5SfaPT+RqMxw0UdHAMY26thLSS2A+/d1DM8D4gzGAxjHM5kHUFTz3knYuQ/7XHcWYGMBvQRh0q4s4NGBhP1HTwaOOW6s4NH/+ZyyvndOTSCgMKx0mUCFgOvGwyGjtD0mqLR96AxOlIANNXoVn8oP9Hc+6d3YBTkSuqnHik44jIcLy2IefV0h7EytgPmr019h6eon9eaEBpBe9NUG+cAix3GwdlAwAgpl+es+h40h44UAE29vK3+UH6iqfZhMBjmOK3GHWEEdBk50wGTixaywbF/lYtlPVrlEm1NU9/hKpf9HeVI1uRzduL4LjvE1d2hHRk8tCTnc27oPWiSDvUDcIxMubgsrxgMhh1Go3FsQ/sDpX2OL/tDxLwwlvqwaIkLzXzGJcB5HaVJNaON8xz7YztqGbCtkI5AEkkQI42AEkkQIwWARBLESAEgkQQxUgBIJEGMFAASSRAjBYBEEsRIASCRBDH/H3iOnduYB2s2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Put model & likelihood into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Initalize plot\n",
    "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "# The gpytorch.settings.fast_pred_var flag activates LOVE (for fast variances)\n",
    "# See https://arxiv.org/abs/1803.06058\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51).cuda()\n",
    "    prediction = likelihood(model(test_x))\n",
    "    mean = prediction.mean\n",
    "    # Get lower and upper predictive bounds\n",
    "    lower, upper = prediction.confidence_region()\n",
    "\n",
    "# Plot the training data as black stars\n",
    "ax.plot(train_x.detach().cpu().numpy(), train_y.detach().cpu().numpy(), 'k*')\n",
    "# Plot predictive means as blue line\n",
    "ax.plot(test_x.detach().cpu().numpy(), mean.detach().cpu().numpy(), 'b')\n",
    "# Plot confidence bounds as lightly shaded region\n",
    "ax.fill_between(test_x.detach().cpu().numpy(), lower.detach().cpu().numpy(), upper.detach().cpu().numpy(), alpha=0.5)\n",
    "ax.set_ylim([-3, 3])\n",
    "ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
